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Abstract 

The issue of developing effective and robust schemes to implement 
a class of the Ogden- type hyperelastic constitutive models is addressed. 
To this end, special purpose functions (running under MACSYMA) are 
developed for the symbolic derivation, evaluation, and automatic FOR- 
TRAN code generation of explicit expressions for the corresponding stress 
function and material tangent stiffness tensors. These explicit forms are 
valid over the entire deformation range, since the singularities resulting 
from repeated principal-stretch values have been theorectically removed. 
The required computational algorithms are outlined, and the resulting 
FORTRAN computer code is presented. 


1 Introduction 

To a great extent, constitutive models of the so-called generalized Rivlin-Mooney 
type [1,2] (i.e., with the stored strain energy density written as a polynomial 
function in terms of the deformation invariants) have dominated the phenomeno- 
logical theory of isotropic hyperelasticity [1-6]. Such models dominate the re- 
lated computational literature on finite-strain elasticity [7-9] as well. Recently 
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though, alternative representations in terms of the principal stretches have be- 
come increasingly popular in nonlinear finite element analyses [6,8,10]. How- 
ever, from the viewpoint of numerical implementation, the use of these models 
presents a number of unique and difficult problems, which do not arise in clas- 
sical representations using the strain invariants. The main difficulty is that (in 
addition to being reasonably complicated functions of the strain components) 
taken separately, the main constituents of the deformation tensor (i.e., prin- 
cipal values and associated eigenvectors) are, in general, not uniquely defined 
and continuously differentiable functions. A careful consideration is thus called 
for in implementing constitutive models formulated in terms of these principal- 
strain measures; this was the main problem addressed by Saleeb and Arnold 
[11] . They bypassed the difficulty entirely by resorting to explicit derivations 
of appropriate forms of the material tangent-stiffness matrices, which are valid 
for the entire deformation range. The explicit expressions they developed [11] 
were for two specific forms of the Ogden-type, strain-energy functions, which 
actually encompass many of the popular representations currently in use for rub- 
ber materials. Results were obtained by simply applying a systematic limiting 
procedure for one type of tensoT- valued function and its spectral representation. 

Symbolic computation specializes in exact computation with numbers, for- 
mulas, vectors, matrices, equations and the like. Numerical computation, on 
the other hand, uses floating-point numbers to compute approximate solutions 
to problems of practical interest. The two approaches are complementary and, 
when combined into an integrated form, can be very powerful in engineering 
applications. In particular, application of symbolic manipulation can provide 
significant incentive for the development of new constitutive theories and their 
applications, for example, finite element. Recently, a problem-oriented, self- 
contained, symbolic expert system, named SDICE (see [12-13]), was developed; 
it is capable of efficiently deriving, in analytical form, potential based constitu- 
tive models whose representations are in terms of the classical invariant formu- 
lation [14-15]. In addition, the FORTRAN code associated with the resulting 
analytical expressions can be automatically generated. 

The objective of the present paper is to discuss three special purpose func- 
tions (SDIFF, SDIFFEV, and TEMPLATE) running under DOE MACSYMA 
[16]. These three functions have been developed to allow the derivation and 
automatic FORTRAN code generation of alternative potential based constitu- 
tive models composed of principal values and their associated eigenvectors, as 
discussed in reference 11. All three functions are written at the MACSYMA 
command level. In the future, these functions will be integrated into the col- 
lection of special purpose functions known as SDICE. This paper begins by 
reviewing highlights of the previous theoretical development and discussing the 
associated computer algorithm for the derivation of the explicit expressions for 
the second Piola Kirchhoff stress tensor Sij and the material moduli tensor 
Dijki • The paper concludes with the evaluation of a separable strain energy 
function, similar to that discussed in reference 11, and its associated FORTRAN 
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source code generation. 


2 Background 

The theoretical development of singularity-free representations for principal 
value-based constitutive models has been discussed at length in reference 11 . 
For brevity, we will confine our discussion, for illustrative purposes, to hypere- 
lastic isotropic materials whose strain energy function W can be taken to have 
the following separable functional dependence: 

W = W(Xi) = J2 a n (A?* + A?* + A 3 *) (1) 

n= 1 

where A,* represents the principal values of the right Cauchy-Green deformation 
tensor C t; . Denoting n t ( i = 1, 2, or 3) to be the associated eigenvectors of 
we can define 

C« = i><i)Af (2) 

7=1 

where is defined as 

( 3 ) 

and is often referred to as the (orthogonal) eigenprojection operator related to 
the associated eigenvectors of Co- 

Equation (2) is valid for the case when all three eigenvalues, A,, are distinct. 
However, for the case when two eigenvalues are the same (i.e., double coalescence 
Ai ^ A 2 = A 3 = A) we have 

Cij = (Ai - A )N§' ) + A 6ij (4) 

And for the case of triple coalescence (Ai = A 2 = A 3 = A), we have 

C i; = A^^' ) =A5 ti . (5) 

7=1 

Similarly, by manipulating equations (2) and (4), we can obtain explicit 
expressions for N$p in terms of Ciji 

- (A r -A.)(A.-A.) ^ * X - S ‘* C ‘> - *** (6) 

and 
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( 7 ) 


N (r) = — ACa - A Sij) 

'* (A r -A) v t; lJ) 

In the preceding equations, the r, s, and t are any cyclic permutation of (1,2, 
or 3). These definitions, equations (2) through (7), will prove very useful in 
obtaining the pertinent singularity-free directional derivatives of both the strain- 
energy potential function W and the stress function Sij = Sij{Cij ). 

The explicit singularity-free expressions for the second Piola Kirchhoff stress 
tensor Sij(Cij ) are defined as 

Sij =2^~ = S ij (C ij ) ( 8 ) 

and those for the material moduli tensor Dijki{Cij ) are obtained by applying 
the directional derivative formula to Sij, that is 


A,h = 2 ^ = 4 


dCijdCti 


= Dijki(Cij) 


As a result, the explicit expressions of the functional dependence of tensors Sij 
and Dijki on Cij can be obtained directly for the following three cases: 1 ) all 
three eigenvalues are distinct; 2) a single singularity (Ai ^ A 2 = A 3 = A, i.e., 
double coalescence) is present; or 3) a double singularity (A, ^ A 2 = A 3 = A, 
i.e., triple coalescence) is present. 


3 Computer Algorithm 

The objective of the present study was to construct three special purpose func- 
tions (SDIFF, SDIFFEV, and TEMPLATE) written at the MACSYMA com- 
mand level that can, respectively, 

(1) Derive explicit expressions for the stress tensor 5», (eqs. ( 8 )) and material 
tensor Dijki (eqs. ( 9 )) ) given three, one, or no distinct eigenvalues 

(2) Evaluate symbolically the expressions generated by SDIFF for a given 
strain-energy function W 

( 3 ) Evaluate the expressions generated by SDIFF and use the built-in MAC- 
SYMA function gentran to automatically generate the associated FOR- 
TRAN code needed to evaluate the expressions numerically for a given 
function, W. 

These special purpose functions contain a list of built-in MACSYMA in- 
structions (factor, expand, ev, ratsubst, diff, limit, and for-loops, to 
name a few) arranged in a specific algorithmic order. Each function, then, can 
be thought of as a macro command. 
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3.1 SDIFF(case) 

Issuing the command SDIFF invokes the following algorithm (consisting of 15 
steps) for automatic derivation of Sij and Dijki ■ In this context, case = 1 
indicates that all three eigenvalues are distinct; case = 2 indicates that only one 
is distinct; and case = 3, that none are distinct. 

To obtain Sii, 

(1) Differentiate W with respect to C,j (see eq. (8)) 


c -V'o 

j 'k dx w 9C <> 


( 10 ) 


(2) Apply the special directional derivative rules obtained from equation (2), 
that is, 

(ii) 


CO _ d\ix 

ij ~ dCa 


whose value is given in equation (6). 

(3) Obtain typical scalar derivatives by using the built-in diff command: 


<‘ 2 > 

(4) Multiply the results s( A(/)) and then sum and factor out coefficients 
of like terms (i.e., CikCkj , Cij, and y ) , thereby obtaining the functional 
dependence of Sij on Qj. In the case of three distinct eigenvalues, 


Sij — aCikCkj -b bCij + cSfj (13) 

where 6ij is the second order identity tensor and 

a = — m[s(A!)(A 2 — A 3 ) 4- s(A 2 )(A 3 — Ai) 4- s(A 3 )(Ai — A 2 )] (14) 


c = 
and 


b = m[s(Ai)(A2 — A3) + $(A 2 )(A3 — Af) + $(^3)(Ai A2)] (15) 

— m[s(A!)A2A 3 (A2 — A 3 ) + s(A 2 )A 3 Ai(A 3 — Ai) 4- $(A 3 )AiA 2 (Ai — A 2 )] 

(16) 


1 

" (Ai - A 2 )(A 2 - A 3 )(A 3 - Ai) 


(17) 


5 



To obtain Dijki • 

(5) Differentiate Sij with respect to Cki (see eq. (9)): 


Dijki = + Sii6 m k)C m j + - Ci m (6jk$mi + 

+ E ^ WT C imC mj + ^(^ + Mj*)l 


r=a 


3 56 5A r ^ ^ dc dX r £ ^ 

°0 + ~ TrT. 0 ^ * 


+ E 5A r dC t r' 3 T ^ 5A r ac if 

( 6 ) Apply the special directional derivative rule 


(0 _ d*(0 


nv> = 

' 3 dC, 


•} 


(7) Obtain the nine scalar derivatives, 


da db dc 

dK’d^’dK 


(18) 


(19) 


(20) 


of equations (14) to (16) for r = 1,2, and 3. 

( 8 ) Substitute the preceding expressions and group-like terms, thus giving 


Dijti = 2aiCkiC?j + 2 a 2 (Ci;Cf ; + C^Cij) + 203 ( 6*1 C?- + C kI 6ij ) 
+2o4(C'fc|Cy) + 2as(6{kCij •+■ CikSj, + SkiCtj -h Gt i Sij ) 

2ae(Sikbji + bki&ij) ( 21 ) 

(9) For comparison of equation (21) to the forms described in reference 11, 
section 4 , we make use of the symmetry properties of Cij and Sij , and 
define two second order symmetric tensors, P and Q, 

PijUG, H) = GaHj, + GuHjk (22) 

Qij kl (G, H) = GaHj, + Gij Hjt + GjiH ik + G jk Hi, (23) 

such that upon substitution we obtain 


Di jkt = a x P(Cl, Cjj) + a 2 [P(C 2 kl , C 0 ) + P(C kh C? )] 
+a 3 [Q(ClSij) + P(6 kl ,Cfj)] + a 4 P(C kl , Cij) 


6 



(24) 


+ &s[Q(Ckh $ij) + Q(&kh C»j)] + 2a$Iijki 

where al, a2, ...a6 are as defined in reference 11 and the preceding equation 
(eq. ( 24 )) is directly comparable to equations 4 . 6 a in reference 11 . Note 
that 


Iijki = + Mi*] ( 25 ) 

Cl = Ci m C m j ( 26 ) 

in the foregoing expressions. 

Next, given the case of nondistinct eigenvalues, for example, case II when (Ai ^ 
A2 = A3 = A), or case III when (Ai = A2 = A3), we must 

( 10 ) Remove the singularity (case II) or singularities (case III) by defining an 
appropriate ”path” for taking the limit of a,b,c and CikCkj in equations 
( 13 ); that is, 

• For case II 

Ai, A2 = A + A, A3 — A — A 

• For case III 

Ai — A, A2 — A + A, A3 = A — A 

(11) Substitute the preceding eigenvalues into the expressions for a,b,and c in 
equations ( 13 ), and take the limit of the numerator and denominator of 
a,b, and c as A— ► 0 . 

(12) If both limits are zero, apply l’HospitaTs rule recursively to the now equiv- 
alent one dimensional problem. For example, given case II, we obtain 


o<*(A) = 




lim^oc(A) = 


1 


JT 1 - KM- s(X) -(Ax -Ay (A)] (27) 

(Aj - A) 2 

[— 2A[s(Aj) - *(A)] + (A? - A 2 )s'(A)] (28) 


(Ai — A ) 2 


1 


(Ai-A) 2 


[A 2 s(Ai) + Ai(Aj - 2A)s(A) - AiA(A x - A)s'(A)] 

(29) 


where 


,,, x _ Hhil - 2W 

1 rJ " 3A (r) SA (r) dA (r ) 


(30) 
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( 31 ) 


(13) Simplify CihCkj by using the definition of Cij and Nij , that is, 

Cf j = X 1 N^ ) + X 2 N^ + X 3 N^ 

In addition, by using Sij = A r // ) + Nfp + N$p and equation (4), for case 
II we obtain 

CtkCki = - X 7 )Cij + (A 2 A a - XXf)Sij] (32) 

and with equation (5) for case III we have 


CitCtj = X6ij. (33) 

(14) Substitute the limiting values of a,b,c and CitCtj into equations (13) and 
group like terms to obtain the modified stress function, Sij , and the a and 
b values for case II: 


where 


and 


For case III, 


Sij — oCij + bS{j 


_ s(Ai) - s(A) 

°“ (A, -A) 


t [«(Ai)A - «(A)Ax] 
(A: -A) 


Sij — s(X)Sij 


(34) 

(35) 

(36) 

(37) 


(15) Repeat steps 5 through 10, but now use the appropriate modified stress 
function. For case II, this results in, 


-«&•*$>+ 1 )<?« + «»» + 

1 (38) 

And for case III, 


da d\ p . 

Dijkl ~ 2 d\ p dCti 6ij 

where the special derivative rule of equation (7) is now used. 


(39) 
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The value in automating the foregoing procedure is evident: not only does 
this special purpose function relieve the user of the tedious manual derivation 
process but it also ensures analytical accuracy. This was illustrated prior to the 
publication of reference 11 in that a number of errors in the hand derivation were 
detected, verified and corrected. Furthermore, as will be discussed in a sequel 
paper [17], this automated derivation procedure facilitated the generalization of 
the preceding expressions to the general nonseparable case, which to the author’s 
knowledge, has eluded researchers to date . Also, it should be apparent that 
this derivation process needs to be executed only once. However, with each 
new definition of W evaluation of s(A(/)) and $'(A(j)) is required in order to 
specialize the needed coefficients; for example, a,b and c, and al, a2, ...a6. As a 
consequence, this motivated the development of SDIFFEV, as described in the 
next section. 


3.2 SDIFFEV (case, W) 

The function SDIFFEV symbolically evaluates the explicit expressions for the 
stress function Sij and material moduli tensor Dijkh which were generated by 
SDIFF and stored in a LISP [18] level disk file. Only the coefficients of these 
expressions need be changed when a different strain-energy function is specified. 
The evaluation algorithm is illustrated here in pseudo code: 

SDIFFEV (case, W) 

IF (diff(W,A 1 ,A 2 ),diff(W,A 2 ,A 3 ), diff(W,A 3} A:))=0 THEN 
Display message: W is separable. 

SEP = 1 

ELSE Display message: W is non separable. SEP = 2 
ENDIF 

IF case=l THEN . 

Call Subroutine A 
ELSE IF case = 2 THEN 
Call Subroutine B 
ELSE IF case = 3 THEN 
Call Subroutine C 
END IF 
End 

Subroutine A 
IF SEP = 2 THEN 
Do loop i — 1, 6 

a[i] = ea[i] (ea[i] are the coefficients of tensor D stored on the disk file produced 
by SDIFF(l)) 

End loop 

ELSE IF SEP = 1 THEN 
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s[2,l] = s[3,l] = s[3,2] = 0 
ENDIF 

Do loop i = 1 , 6 
a[i] = ev(ea[i]) 

End loop 
Do loop i = 1 , 3 
s[i] = 2*diff(W,Ai,l) 
s[i,i] = 2*difF(W,Aj,2) 

IF SEP = 2 THEN 
Do loop j = 1, 3 
s[ij] = diff(W,A,, A ; ,2) 

End loop 
ENDIF 
End loop 
Call OPTION 
Return End 

Subroutine B 
W = ev(W,A 3 = A 2 ) 

IF SEP = 2 THEN 
Do loop i = 1, 3 

b[i] = eb[i] (eb[i] are the coefficients of tensor D stored on the disk file produced 
by SDIFF (2)) 

End loop 

ELSE IF SEP = 1 THEN 

s[2,I] = 0 

Do loop i = 1, 6 
b[i] = ev(eb[i]) 

End loop 
Do loop i = 1, 2 
s[i] = 2*diff(W,Ai,l) 
s[i,i] = 2*diff(W,A„2) 

IF SEP = 2 THEN 
Do loop j = 1, 2 
s[ij] = diff(W,A,, A ; ,2) 

End loop 
ENDIF 
End loop 
Call OPTION 
Return 
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Subroutine C 
W = ev(W,A 3 = A 2 = Ai) 
s[l]=2*diff(W,A 1 ,l) 
s[l,l]=2*diff(W,A 1 ,2) 

Call OPTION 
Return 

Subroutine OPTION 

Display the formulae S[ij] and D[ij,k,l]. Then, ask if user wants to see the 
symbolic form for the given function W, the intermediate step evaluations, and 
the derivatives of W. 

READ(type y, or n to the question) 

DISPLAY the options user may choose 
Return 


3.3 TEMPLATE () 

The function TEMPLATE is similar to the function SDIFFEV in that both 
will evaluate the explicit expressions obtained from SDIFF . As a result neither 
can be employed unless preceded by an invocation of SDIFF. TEMPLATE, 
however, will automatically generate the associated FORTRAN source code 
needed to evaluate the expressions numerically for a given potential function W. 
Code generation is accomplished by utilizing the built-in MACSYMA function 
gentran, and a number of template files. The template files can be thought of 
as a framework for the generation of four basic FORTRAN subroutines (i.e., the 
main driving routine COMPSD and the three subroutines — one each for case I , 
case II, and case III) and numerous functions. Appendix A contains the template 
file for the main driving routine COMPSD. This subroutine is constructed for 
easy implementation into a finite element code; the input requirements are the 
strain tensor C m (denoted as emu) and its associated eigenvalues (i.e., Ai j A 2 , A 3 
denoted by gll, gl2, and gl3 respectively), and the outputs are the stress tensor 
S n (denoted as s), and the material moduli tensor D nm (denoted as d). Here, 
71 and m run from 1 to 6. The only automated code generation required is that 
for the subroutines COMPSD1, COMPSD2, and COMPSD3. These codes are 
generated by issuing the command <gentranin> . The subroutines COMPSD1, 
COMPSD2, and COMPSD3 are associated with case I (Ai ^ A 2 ^ A 3 ), case II 
(Ai, A = A 2 = A 3 ), and case III ( A = X x = A 2 = A 3 ), described in Section 2.0. 
The template files corresponding to these three cases are shown, respectively in 
appendixes B,C, and D. Note that in these routines, most of the FORTRAN 
code is automatically generated, since it pertains to the definition of coefficients 
a,b,c; al, a2, ..., a6, and the first (sl,s2,$3, see eq. (12)) and second scalar 
(sll,s22,s33 , see eq. (30)) derivatives of the strain energy function W. The 
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gent ran commands are enclosed by double inequality signs, that is, ^ > * 
Finally, all functions that are associated with a given case have been included 
in the corresponding appendix. 

4 Example 

As an example, consider the case in which the strain energy function W of 
equation (1) consists of only two terms; that is, 

W = xl(gll yl + gl 2 yl + gl 3 yl ) + x2(gll y2 + gl2 y2 + glZ y2 ) (40) 

where xl, x2, yl, and y2 are material coefficients and gl 1 = A li gl2 = A 2 , 
and glS = A 3 . After defining W, we can symbolically obtain the analytical 
expressions for Sij and D,j ki (given the case of three distinct eigenvalues) by 
merely issuing the command 


sdiffev(l, W)\ 

at the MACSYMA command level. Case II or III can just as easily be obtained 
by substituting a 2 or 3 in place of the 1 in this command. The resulting output 
is shown in appendix E where the expressions for the coefficients a,b,c and 
al ,a2,. . .a6 could be further simplified and manipulated, if desired, by using other 
MACSYMA built-in functions. Typically, however, the analyst will ultimately 
desire a FORTRAN code for the resulting expressions in order to solve a given 
structural problem using the foregoing constitutive model. This code, described 
in the previous section, can easily be obtained by issuing the command 

template(); 

at the MACSYMA command level. The generated FORTRAN code will then 
be stored in a file named temp.f. The automatically generated FORTRAN code 
for the above example is shown in appendix F. 

5 Summary of Results 

Taken separately, the main constituents of the deformation tensor (i.e., prin- 
cipal values and associated eigenvectors) are, in general, not uniquely defined 
and continuously differentiable functions. Careful consideration is thus called 
for in implementing constitutive models formulated in terms of these principal- 
strain measures. This difficulty can be entirely bypassed by resorting to ex- 
plicit symbolic derivations of appropriate forms of the material tangent-stiffness 
matrices, which are valid over the entire deformation range. Furthermore, to 
enhance effective utilization and implementation of the present results, auto- 
matic FORTRAN generation of these explicit expressions has been pursued and 
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achieved. As a result, three special purpose functions(SDIFF, SDIFFEV and 
TEMPLATE), running under MACSYMA, have been developed and verified. 


References 

[1] Green, A.E. and Adkins, J.E.: Large Elastic Deformations . 2nd Edition. 
Oxford TJniv. Press, Clarendon, 1970. 

[2] Rivlin, R.S.: Large Elastic Deformations of Isotropic Materials; Part III: 
Experiments on the Deformation of Rubber. Philo. Trans. Roy. Soc. , Lon- 
don, Ser. A, Vol. 243, pp.251-288, 1951. 

[3] Rivlin, R.S.: ’’Forty Years of Nonlinear Continuum Mechanics”, 

Proc. IX Int. . Congress on Rheology, Mexico, 1984. 

[4] Treloar, L.R.G.: The Physics of Rubber Elasticity , 3rd Ed., Oxford Univ. 
Press., U.K., 1975. 

[5] Ogden, R.W., ’’Recent Advances in the Phenomenological Theory of Rub- 
ber Elasticity,” J.Rubber Chem. Tech- Vol. 59, pp. 361-383, 1986. 

[6] Finney, R.H., and Kumar, A.: "Development of Material Constants for 
Non-linear Finite Element Analysis,” J. Rubber Chem. Tech. . Vol.. 61, pp. 
879-891, 1988. 

[7] Oden, J.T.: Finite Elements of Nonlinear Continua , McGraw Hill, New 
York, 1972. 

[8] Sussman, T.; and Bathe, K.J.: A Finite Element Formulation for Nonlinear 
Incompressible Elastic and Inelastic Analysis, Comp. Struct. , Vol. 26, pp. 
357-409, 1987. 

[9] Simo, J.C.; Taylor, R.L; and Pister, K.S.: "Variational and Projection 
Methods for the Volume Constraint in Finite Deformation Elastoplastic- 
ity”, Comp. Meth. Appl. Mech. Engng. , Vol. 51, pp. 177-208, 1985. 

[10] Chang, T.Y., Saleeb, A.F., and Li, G.: "Large Strain Analysis of Rubber- 
Like Materials by a Mixed Finite Element”, Computational Mech. , Vol. 8, 
No. 4, pp. 221-233, 1991. 

[11] Saleeb. A. F.; and Arnold, S.M.: Explicit Robust Schemes for Implementa- 
tion of a Class Of Principle Value-Based Constitutive Models: Theoretical 
Development. NASA TM 105345, 1991 

[12] Arnold, S.M.; and Tan, H.Q.: "Symbolic Derivation of Potential Based 
Constitutive Equations”, Computational Mech. , Vol. 6, pp. 237-246, 1990. 


13 



[13] Arnold, S.M.; Tan, H.Q.; and Dong, X.: Application of Symbolic Com- 
putations to The Constitutive Modeling of Structural Materials. Symbolic 
Computations and Their Impact on Mechanics , Noor, A.K., Elishakoff, I. 
and Hulbert, G., eds., PVP-Vol. 205, ASME, pp. 215-229. ,1990. 

[14] Malvern, L.E.: Introduction to the Mechanics of_a Continuou s Medium. 

Prentice- Hall, 1969. 

[15] Chen, W.F.; and Saleeb, A.F.: Constitutive Equations For Engineering 
Materials Volume 1: Elasticity and Modeling , John Wiley & Sons, 1982. 

[16] MATHLAB Group: MACSYMA Reference Manual . Version 10. Labora- 
tory for Computer Science, Massachusetts Institute of Technology, Cam- 
bridge, MA, 1984. 

[17] Arnold, S.M.; et al.: Explicit Robust Schemes for Implementation 

of General Principal Value-Based Constitutive Models, submitted to 
Comp, and Struct. , 1993. 

[18] Wilensky, R.: LISPcraft . W.W. Norton Co, 1984. 


14 



APPENDIX A: Template File Associated With COMPSD 
The Main Driver Routine 

c **************** ******************** ***************** 
c This is the template subroutine to calculate 

c tensor S and D. inputs are eigenvalues gll,gl2,gl3, 

c and cmu(6) . emu is assumed to be engineering strain(e), 

c e.g. the Cauchy-green deformation tensor cm(3,3) is related 

c to cmu(6) in the following fashion: 

c cm(l,l)=cmu(l) , cm(2,2) = cmu(2) , cm(3,3) =cmu(3) , 

c cmu(4)=2*cm(l ,2) , cm(5)=2*cm(2,3) ,cmu(6)=2*cm(l ,3) . 

c The outputs are the second order tensor S(6) 

c and forth order tensor D(6,6) are related in the 

c following way: 

c S=D*C 

c S(l,l) = SCD 

c S (2 , 2) = S (2) 

c S(3,3) = S (3) 

c S (1 , 2) = S (4) 

c S (2 ,3) = S ( 5 ) 

c S(3,l) = S(6) 

c C (1 , 1) = C(l) 

c C(2,2) = C(2) 

c C(3,3) = C(3) 

c C(1 ,2) = C(4) 

c C(2 ,3) = C(5) 

c C(3 , 1) = C(6) 

c 

subroutine compsd(gll ,gl2 ,gl3,cmu,s,d) 
real*8 gll,gl2,gl3,ts(3,3) ,td(3,3,3,3) 
real*8 delt(3,3) ,delt4(3,3,3,3) ,s(6) ,d(6,6) 
real*8 cmu(6) ,cm(3,3) 
c 

c converts cmu(6) to matrix cm(3,3) in a way that 

c cm(l,2)=cm(2,l)=cmu(4) ,cm(2,3)=cm(3,2)=cmu(5) , 

c cm(l,3)=cm(3,l)=cmu(6) . 

c 
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do 5 i=l,3 
do 5 j=l,3 

if (i.eq.j) then 
iq=i 

cm(i, j)=cmu(iq) 
else if (i.ne.j) then 
if ((i+j) .eq.3) iq=4 
if( (i+j) .eq.4) iq=6 
if C (i+j) .eq.5) iq=5 
cm(i, j)=cimi(iq)/2 
end if 
continue 

5 continue 
c 

c Initiates the second identity tensor delt(3,3) which 

c is a 2X2 identity matrix, 

c 

do 6 i=l,3 

delt(i,i)=1.0 

6 continue 
c 

c Computes the forth order identity tensor delt4(3,3,3,3) 

c by definition, 

c 

do 7 i-1,3 
do 7 j=i,3 

delt4(i , j ,i, j)=delt(i,i)*delt(j , j)+delt(i, j)*delt(j ,i) 
delt4(i, j , j ,i)=delt4(i, j ,i, j) 

7 continue 
c 

***************** ****************************************** ******* 
c For different eigenvalues gll,gl2,gl3 the computation 

c is different, easel is gll#gl2#gl3 call subroutine comsdl . 

c case2 is gl3=gl2#gll or gll=gl3#gl2 or gll=gl2#gl3 then 

c call subroutine compsd2 . case3 is gll=gl2=gl3 call subroutine 

c compsd3 . 

c ******************************************************************** 
if ((gll.ne.gl2) .and. (gl2.ne.gl3) .and. (gll.ne.gl3)) then 
call compsdl(gll,gl2,gl3,delt ,delt4,cm,ts > td) 
else if (Cgl2.eq.gl3) .and. (gll.ne.gl3)) then 
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call compsd2(gll,gl2,delt,delt4,cm,ts,td) 
else if ((gll .eq.gl2) .and. (gl3.ne.gl2)) then 
gll=gl3 

call compsd2 (gll , gl2 , delt , delt4 , cm, ts , td) 
else if ( (gll . eq.gl3) .and. (gl2 .ne .gl3)) then 
gll=gl2 
gl2=gl3 

call coinpsd2(gll ,gl2,delt ,delt4,cm,ts ,td) 

else 

call coinpsd3 (gll , delt ,delt4,ts ,td) 
end if 
c 

c Rewrite the tensor ts(i,j) td(i, j ,k,l)to S(i) and D(i,j) 

c respectively by using the symetric property, 

c converts ts(3,3) s(6) and td(3,3,3,3) to D(6,6) 

c 

do 8 i=l,3 
do 8 j=i,3 

if (i.eq.j) iq=i 
if (i.eq.l.and. j .eq.2) iq=4 
if (i. eq.2. and. j .eq.3) iq=5 
if (i . eq. 1 .and. j . eq.3) iq=6 
s(iq)=ts(i, j) 
cont inue 

8 continue 
do 9 i=l,3 

do 9 j=i,3 

d(i , j)=td(i , i , j , j) 
continue 

9 continue 

do 10 i=l ,3 

d(i ,4)=td(i,i , 1 ,2)+td(i ,i,2 ,1) 
d(i , 5)=td(i ,i ,2 ,3)+td(i , i,3 ,2) 
d(i,6)=td(i,i,3, l)+td(i ,i, 1 ,3) 

10 continue 

d(4,4)=(td(l , 2 , 1 , 2)+td(l , 2 , 2, l)+td(2, 1 ,1 ,2) +td(2 , 1 ,2 , l))/2 
d(4, 5)=(td(l ,2,2,3)+td(l,2,3, 2)+td(2 ,1,2,3) +td(2 , 1 ,3, 2) ) /2 
d(4 ,6)=(td(l ,2,l,3)+td(l > 2,3, l)+td(2 ,1,1,3) +td(2 , 1 ,3, l))/2 
d(5,5)=(td(2,3,2,3)+td(2,3 > 3,2)+td(3 > 2,2,3)+td(3 > 2,3,2))/2 
d(5 ,6)=(td(2 ,3, 1 ,3)+td(2,3 ,3, l)+td(3,2 ,1 ,3)+td(3,2,3, l))/2 
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d(6 ,6)=(td(3, 1 , 1 ,3)+td(3 , 1 , 3 , l)+td(l,3,l ,3)+td(l ,3,3, l))/2 . 
do 11 i = 1,6 
do 11 j = 1,6 

d(i , j) = d(j , i) 

11 continue 

c 

c prints out the inputs gll ,gl2,gl3,cniu(6) and outputs S and D 

c 

print*, ’gll 58 ’, gll 

print*, ’gl2=’, gl2 

print*, ’gl3=’, gl3 

print*, ’Input tensor C(6):’ 

print*, (cmu(i), i * 1,6) 

print*, "second order tensor S(6):" 

print*, (s(i), i=l,6) 

print*, "The forth order tensor D(6,6):" 
do 101 i=l ,6 

print*, (d(i, j) ,j=l,6) 

101 continue 
return 
end 
c 

subroutine compsdl (gll , gl2 ,gl3 ,delt ,delt4,cm, ts , td) 

« 

gentranin(" easel 1 .tem")$ 

» 

subroutine compsd2 (gll , gl2 , delt , delt4 , cm, t s , td) 

« 

gentranin("case22 . tem") $ 

» 

subroutine compsd3(gll , delt ,delt4,ts,td) 

« 

gentranin("case3 ,tem")$ 

» 
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c 

c This subroutine computes P and Q forth order tensors 

c which we define in tensor D. 

c 

subroutine pqcom(cml,cm2,p,q) 

real*8 cml(3,3) ,cm2(3,3) , p(3,3,3,3) ,q(3,3,3,3) 
do 100 i=l ,3 
do 100 j-1,3 
do 100 k=l ,3 
do 100 1=1,3 

p(i, j ,k,l)=cml(i,k)*cm2(j ,l)+cml(i,l)*cm2(j ,k) 
q(i , j ,k,l)=p(i , j ,k,l)+cml(j ,l)*cm2(i ,k)+cml(j ,k)*cm2(i,l) 
100 continue 

return 
end 
c 

c This subroutine computes matrix product cmXcm. 

c 

subroutine product (matl , cmm) 
real*8 matl (3 ,3) , cmm(3,3) , sum 
do 25 i=l,3 
do 25 j=l ,3 
sum=0 . 0 
do 26 k=l ,3 

sum= sum+mat 1 ( i , k) *mat 1 ( k , j ) 

26 continue 

cmm(i , j)=smn 
25 continue 
return 
end 
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APPENDIX B: Template File Associated With C0MPSD1 
Valid For Three Distinct Eigenvalues 


real*8 gll,gl2,gl3 jts(3,3) ,td(3,3,3,3) 

real*8 cm(3,3) , delt (3,3) , delt 4 (3, 3, 3, 3) ,p(3,3,3,3) 

real *8 q(3 , 3 , 3 , 3) , cinm(3 ,3) , pi (3 ,3 ,3 ,3) ,p21 (3 , 3 ,3 ,3) 

real *8 p31 (3, 3, 3, 3) , qll (3 , 3 ,3 , 3) ,ql2 (3 ,3 ,3 , 3) ,p22 (3 ,3 , 3 ,3) 

real*8 q21(3,3,3,3) ,q22(3,3,3,3) ,a,b,c,ai,a2,a3,a4,a5,a6 

Obtains cinm(3,3) : =cm(3,3)*cin(3,3) from subroutine product 
call product (cm, cmm) 

Uses the formula we derived in code to compute second order 
tensor ts(3,3) . 


gentran(for i:l thru 3 do 
(for j : 1 thru 3 do 

(t s [i , j ] : a (gl 1 , gl2 , gl3) * cmm [i , j ] +b (gll , gl2 , gl3) 
*cm[i, j]+c(gll ,gl2,gl3) *delt [i, j] )))$ 


Call subroutine to compute all the functions we defined 
when we derived forth order tenosor td, namely P(i,j,k,l) 
and Q(i,j,k,l) which are the functions of cm(3,3) and 
the matrix product cmm(3,3). 

call pqcom ( cmm , cmm , pi , q) 
call pqcomCcmm, cm,p21 ,q) 
call pqcom(cm ) cmm,p22 ,q) 
call pqcom(cm,cm,p31 ,q) 
call pqcom(cntm, delt ,p ,qll) 
call pqcom(delt ,cmm,p,ql2) 
call pqcom (cm, delt ,p,q21) 
call pqcom(delt ,cm,p,q22) 



c Computes forth order tensor td(i,j,k,l) 

c 

« 

gentran(for i:l thru 3 do 
(for j : 1 thru 3 do 
(for k:l thru 3 do 
(for 1 : 1 thru 3 do 

(td [i , j , k , 1] : al (gll , gl2 , gl3) *pl [i , j ,k, 1] +a2 (gll , gl2 , gl3) 
* (p21 [i , j , k , 1] +p22 [i , j , k , 1] ) +a4 (gll , gl2 , gl3) *p31 [i , j , k , 1] 
+a3 (gll ,gl2 ,gl3) * (qll [i , j ,k,l] +ql2 [i , j ,k, 1] )+ 
a5 (gll ,gl2 ,gl3)*(q21 [i , j ,k,l] +q22[i,j ,k,l]) + 
a6 (gll , gl2 , gl3)*delt4 [i , j ,k , 1] ) ) ) ) ) $ 

» 

return 

end 

c 

c a,b,c,al-a6 are the coefficients we derived in code, 

c 

« 

gentran(a(gll ,gl2 ,gl3) : =block(type (function, a) , 

typeO'real+S" ,gll,g!2,gl3) , 
type( n real*8 n ,a, si , s2,s3) , 
a:eval(ta)))$ 

» 

« 

gentran(b(gll ,gl2 ,g!3) : =block (type (funct ion, b) , 

type ( n real*8 11 ,b , gll , gl2 , gl3) , 
type( n real*8 n ,sl ,s2,s3) , 
b:eval(tb)))$ 

» 

« 


gentran(c(gll ,gl2,gl3) : =block (type (function, c) , 

type( ,l real*8 lf , c ,gll ,g!2 ,gl3) , 
type( n real*8 M ,sl ,s2 ,s3) , 
c:eval(tc)))$ 
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« 

gentran(al(gll,gl2,gl3) :=block(type(function,al) , 

type("real*8",al,gll,gl2,gl3) , 
type ( "real*8" , si , s2 , s3 , sll , s22 , s33) , 
al : eval(tal)))$ 

» 

« 

gentran (a2 (gll , gl2 , gl3) :=block(type (function, a2) , 
type ("real*8",a2, gll, gl2,gl3) , 
type("real*8" , si ,s2,s3,sll , s22 ,s33) , 
a2 : eval (ta2) ) ) $ 

» 

« 

gentran (a3 (gll , gl2 , gl3) :=block (type (funct ion, a3) , 
type("real*8" ,a3,gll,gl2,gl3) , 
type ("real*8" , si , s2 , s3,sll , s22 ,s33) , 
a3:eval(ta3)))$ 

» 

« 

gentran (a4 (gll , g!2 , gl3) :=block (type (funct ion, a4) , 
type("real*8" ,a4,gll,gl2,gl3) , 
type ("real*8" , si , s2 , s3 , sll , s22 , s33) , 
a4:eval(ta4)))$ 

» 

« 

gentran (a5 (gll ,gl2 ,gl3) :=block (type (funct ion, a5) , 
type ("real*8" , a5 ,gll ,gl2 ,gl3) , 
type ("real*8" , si ,s2,s3,sll , s22 ,s33) , 
a5 :eval(ta5)))$ 

» 

« 

gentran (a6 (gll, gl2,gl3) : =block (type (function, a6) , 
type ("real*8" , a6 ,gll ,gl2 ,gl3) , 
type ( "real*8" , si , s2 , s3 , sll , s22 , s33) , 
a6:eval(ta6)))$ 

» 
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c 


c 

c 


« 


» 


c 


« 


» 


c 


« 


» 


c 


« 


» 


c 


The si , s2, s3, sll , s22, s33 are derivatives of W 

function si (gll ,gl2 ,gl3) 

<<cut (var ) ;» 

gentran(type( "real*8" ,sl,gll,gl2,gl3) , 
si :2*eval(diff (s, 'gll,!)))! 


return 

end 

function s2 (gll ,gl2 ,gl3) 

<<cut (var) ; » 

gentran(type( "real*8" ,s2,gll,gl2,gl3) , 
s2 :2*eval(diff (w, ’gl2,l)))$ 


return 

end 

function s3(gll ,gl2 ,gl3) 

<<cut (var) ; » 

gentran(type( "real*8" ,s3,gll,gl2,gl3) , 

s3 : 2*eval(diff (w, ’gl3,l)))$ 


return 

end 

function sll (gll ,gl2,gl3) 

<<cut (var) ; » 

gentran(type( "real*8" ,sll ,gll,gl2,gl3) , 

sll : 2*eval(diff (w, ’gll ,2)))$ 


return 

end 


23 



function s22(gll ,gl2,gl3) 

<<cut (var) ; >> 

« 

gentran(type("real*8" ,s22,gll ,gl2,gl3) , 
s22 :2*eval(diff (w, ^12,2)))$ 

» 

return 

end 

c 

function s33(gll ,gl2,gl3) 

<<cut (var) ; >> 

« 

gentran(type( "real*8",s33,gll,gl2,gl3) , 
s33 : 2*eval (dif f (w , 1 gl3 , 2) ) ) $ 

» 

return 

end 
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APPENDIX C: Template File Associated With C0MPSD2 
Valid For Double Coalesence Case 


c 

c 

real*8 gll ,gl2,ts(3,3) ,td(3,3,3,3) 
real*8 cm(3 , 3) , delt (3,3), delt4 (3 , 3 , 3 , 3) , pi (3 , 3 , 3 , 3) 
real*8 q2(3,3,3,3) ,ql(3,3,3,3) ,p(3,3,3,3) ,q(3,3,3,3) 
real*8 bl,b2,b3, abar.bbar 
c 

c Computes second order tensor ts(i,j) based on the formula 

c derived in code, 

c 

« 

gentranCfor i:l thru 3 do 
(for j : 1 thru 3 do 

(ts [i , j] :abar(gll ,gl2)*cm[i, j]+bbar(gll ,gl2)*delt [i, j] )) )$ 

» 

c 

c Call subroutine to get P, Q which are defined in code, 

c 

call pqcom(cm,cm,pl ,q) 
call pqcom ( cm , delt , p , ql ) 
call pqcom(delt,cm,p,q2) 
c 

c Computes tensor td(i,j,k,l). 

c 

« 

gentran(for i:l thru 3 do 
(for j : 1 thru 3 do 
(for k : 1 thru 3 do 
(for 1:1 thru 3 do 

(t d [i , j , k , 1] : bl (gll , gl2) *pl [i , j , k , 1] +b2 (gll , gl2) * 

(ql [i , j ,k,l]+q2 [i, j ,k,l] )+b3(gll,gl2)* 
delt4[i, j ,k,l] )))))$ 

» 

return 

end 
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c 

c abar.bbar are bl , b2, b3 functions derived in code, 

c 

« 

gentran(abar(gll ,gl2) :=block (type (function, abar) , 

type("real*8", abar, gll, gl2) , 
type("real*8" , ssl,ss2), 
abar : eval (abar) ) ) $ 

>> 

« 

gentran(bbar(gll ,gl2) :=block(type(function,bbar) , 

type ("real*8",bbar, gll, gl2) , 
type("real*8" , ssl,ss2), 
bbar : eval (bbar) ) ) $ 

» 

« 

gentran(bl (gll ,gl2) :=block (type (functional) , 

type ("real*8",bl, gll, gl2) , 
type("real*8" , ssl,ss2,ssll,ss22) , 
blreval (tbl)))$ 

» 

« 

gentran(b2(gll ,gl2) :=block(type(function,b2) , 

type ("real*8",b2, gll ,gl2) , 
type("real*8" , ssl , ss2 , ssll , ss22) , 
b2 : eval (tb2) ) ) $ 

» 

« 

gentran(b3 (gll ,gl2) :=block(type(function,b3) , 

type("real*8" ,b3,gll ,gl2) , 
type("real*8" , ssl , ss2 , ssll , ss22) , 
b3 : eval (tb3) ) ) $ 

» 

« 

neww : sub st ( [ ’ gl3= 1 gl2] , w) $ 

» 
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c 

c 

c 


« 


» 


c 


« 


» 


c 


« 


» 


c 


<< 


» 


ssl ,ss2 , ssll, ss22 are derivatives of W. 

function ssl(gll,gl2) 

<<cut (var) ;» 

gentran(type("real*8" ,ssl ,gll ,gl2) , 

ssl :2*eval(diff (new, ’gll,l)))$ 


return 

end 

function ss2(gll,gl2) 

«cut (var) ; » 

gentran(type("real*8 M ,ss2,gll,gl2) , 
ss2 :2*eval(diff (neww, J gl2,l)))$ 

return 

end 

function ssll (gll ,gl2) 

<<cut (var) ; » 

gentran(type( M real*8" , ssll , gll ,gl2) , 

ssll :2*eval(diff (new, ’gll, 2)))$ 


retiirn 

end 

function ss22(gll ,gl2) 

«cut(var) ;» 

gentran(type("real*8" , ss22 ,gll ,gl2) , 

ss22 : 2*eval(diff (new, ’gl2,2)))$ 


return 

end 
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APPENDIX D: Template File Associated With C0MPSD3 
Valid For The Triple Coalesence Case 


c 

real*8 gll,ts(3,3) ,td(3,3,3,3) ,delt(3,3) ,delt4(3 ,3,3,3) 
real*8 ccl,abbar 
c 

« 

gentran(for i:l thru 3 do 
(for j : 1 thru 3 do 
(ts[i, j] : abbar (gll )*delt [i, j])))$ 

» 

« 

gentran(for i:i thru 3 do 
(for j : 1 thru 3 do 
(for k:l thru 3 do 
(for 1:1 thru 3 do 

(td[i, j ,k,l] :ccl(gll)*delt4[i, j ,k,l] )))))$ 

» 

return 

end 

« 

gentran(abbar(gll) :=block(type (function, abbar) , 

type("real*8" , abbar, gll), 
abbar : eval (abbar) ) ) $ 

» 

« 

gentran(ccl(gll) :=block(type(function,ccl) , 

type("real*8" , ccl.gll), 
ccl :eval(ccl)))$ 

» 

« 

: subst ( [ } gl3= * gll , ’ gl2= ’ gll] , w) $ 

» 
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c 

function sssl(gll) 

<<cut (var) ; » 

« 

gentran(type("real*8" ,sssl ,gll) , 
sssl : 2*eval(diff (www, ^11,1)))$ 

» 

return 

end 

c 

function sssll(gll) 

<<cut (var) ; » 

« 

gentran(type("real*8" ,sssll ,gll) , 

sssll :2*eval(diff (www, *gll,2)))$ 

» 

return 

end 


29 



APPENDIX E: Listing of MACSYMA Session Resulting From 
Issuing The SDIFFEV Command 


c6) sdiffev(l,w) ; 

w is a separable function. 

y2 y2 y2 yl yl yl 

w = (gl3 + gl2 + gll ) x2 + (gl3 + gl2 + gll ) xl 

This is case 1 with distinct eigenvalues gll#gl2#gl3. 

Please type y if your answer is yes, otherwise type n to skip it. 

Do you want to display the second order tensor s[i,j]? 

V> 

s - a c c + c delt + b c 
i, j i, k k, j i, j i, j 

Do you want to display c[i,j] and delta[i,j]? 

y; 

[10 0] 

[ ] 

delt =[010] 

i » j [ 3 

[ 001 ] 

c = gl3 n3 n3 + gl2 n2 n2 + gll nl nl 

i> j i j i j i j 

nl, n2, n3 axe eigenvectors associatedd with eigenvalues gll, gl2, gl3. 
If c[i,j]) is given then the eigenvectors can be computed 
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Do you want to display a,b,c in s[i,j] form? 

y; 

s3 s2 si 

a = + 

t31 (t31 - t21) t21 (t31 - t21) t21 t31 

s3 (2 t31 - t21 -2 gl3) si (t31 + t21 + 2 gll) 

b = + 

t31 (t31 - t21) t21 t31 

s2 (t31 - t21 + gl2 + gll) 
t21 (t31 - t21) 

s3 (t31 - gl3) (t31 - t21 - gl3) s2 (t21 - gl2) (t31 + gll) 

c = 

t31 (t31 - t21) t21 (t31 - t21) 

si (t21 + gll) (t31 + gll) 


t21 t31 

Do you want to display t21,t31, sl,s2,s3? 

y; 

t21 = gl2 - gll 
t31 = gl3 - gll 

sl,s2,s3 are the first derivatives of W with respect to gll,gl2,gl3. 

y2 - 1 yl - 1 

si = gll x2 y2 + gll xl yl 
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y 2 - 

1 

yl - 

1 

s2 = gl2 

x2 y2 + gl2 

xl yl 

y2 - 

1 

yl - 

1 

s3 = gl3 

x2 y2 + gl3 

xl yl 


Do you want to display the forth order tensor d[i,j,k,l]? 

y; 


2 2 

d * [a delt4 + a (q(delt, c ) + q(c , delt)) 


+ a (q(delt , c) + q(c, delt)) + a p(c , c ) + a (p(c , c) + p (c, c )) 
5 1 2 

+ a p(c, c)] 

4 

Do you want to display the functions p,q and delt4? 

y; 


p(g, h) = g h +g h 

i , k j,l i . 1 j i ^ 


q(g> h)=g h +h g +2g h 

i, k j, 1 i, k j, 1 i, 1 j, k 


delt4 = delt delt + delt delt 

i,k j , 1 i.l j>k 


Do you want to continue displaying a 


y; 
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2 s3 (2 t31 - t21) 2 si (t31 + t21) s33 . 


1 3 3 3 3 2 2 

t31 (t31 - t21) t21 t31 t31 (t31 - t21) 

s22 2 s2 (t31 - 2 t21) sll 


2 2 3 3 2 

t21 (t31 - t21) t21 (t31 - t21) t21 


Do you want to continue displaying a ? 

2 

y; 

2 2 
si (2 t31 + 3 t21 t31 + 4 gll t31 + 2 t21 + 4 gll t21) 


2 3 3 

t21 t31 

2 2 
s2 (2 t31 - 3 t21 t31 + 4 gll t31 - t21 - 8 gll t21) 


3 3 

t21 (t31 - t2l) 

2 2 
s3 Ct31 + 3 t21 t31 + 8 gll t31 - 2 t21 - 4 gll t21) 


3 3 

t31 (t31 - t2l) 

s33 (2 t31 - t21 - 2 gl3) sll (t31 + t21 + 2 gll) 


2 2 2 2 
t31 (t31 - t2l) t21 t31 


2 

t31 
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s22 (t31 - t21 + gl2 + gll) 


2 2 
t21 (t31 - t21) 

Do you want to continue displaying a ? 

3 

y; 

3 2 2 2 2 
a - - s3 (t3i - 2 t21 t31 - gll t31 + t21 t31 - 3 gll t21 t31 -4 gll t31 

3 

2 2 3 3 

+ 2 gll t21 + 2 gll t21) / (t31 (t31 - t2l) ) 

2 2 2 2 

- si (t21 t31 + 2 gll t31 + t21 t31 + 3 gll t21 t31 + 2 gll t31 

2 2 3 3 

+ 2 gll t21 + 2 gll t21)/(t21 t31 ) 

2 2 2 2 
+ s2 (t21 t31 + 2 gll t31 - 2 t21 t31 - 3 gll t21 t31 + 2 gll t31 + t21 

2 2 3 3 3 

- gll t21 - 4 gll t2l) / (t21 (t31 - t21) ) 

s33 (t31 - gl3) (t31 - t21 - gl3) s22 (t21 - gl2) (t31 + gll) 


2 2 2 2 
t31 (t31 - t2l) t21 (t31 - t2l) 

sll (t21 + gll) (t31 + gll) 


2 2 
t21 t31 
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Do you want to continue displaying a ? 

4 

y; 

2 2 

2 s3 (t21 + 2 gll) Ct31 + t21 t31 + 4 gll t31 - t21 - 2 gll t21) 

a 

4 3 3 

t31 (t31 - t21) 

2 2 

2 si (t31 + t21 + 2 gll) (t31 + t21 t31 + 2 gll t31 + t21 +2 gll t21) 


3 3 

t21 t31 

2 2 

2 s2 (t31 + 2 gll) (t31 - t21 t31 + 2 gll t31 - t21 - 4 gll t21) 



3 3 

t21 (t31 - t2l) 

2 2 
s33 (2 t31 - t21 - 2 gl3) sll (t31 + t21 + 2 gll) 


2 2 2 2 
t31 (t31 - t21) t21 t31 

2 

s22 (t31 - t21 + gl2 + gll) 


2 2 
t21 (t31 - t21) 


Do you want to continue displaying a ? 

5 


y; 
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2 


3 3 2 2 2 2 

a = si (t21 t31 + 2 gll t31 + t21 t31 + 6 gll t21 t31 + 6 gll t31 

5 

3 2 2 3 3 

+ t21 t31 + 6 gll t21 t31 + 9 gll t21 t31 + 4 gll t31 + 2 gll t21 

2 2 3 3 3 

+ 6 gll t2l + 4 gll t2l)/(t21 t31 ) 

3 3 2 2 2 2 2 

+ s3 (t21 t31 + 2 gll t31 - 2 t21 t31 - 6 gll t21 t31 - 3 gll t31 

3 2.3 3 2 2 

+ t21 t31 - 9 gll t21 t31 - 8 gll t31 + 2 gll t21 + 6 gll t21 

3 3 3 

+ 4 gll t2l) / (t31 (t31 - t21) ) - s2 

3 3 2 2 2 2 3 

(t21 t31 + 2 gll t31 - 2 t21 t3i + 6 gll t31 + t21 t31 - 

2 2 3 3 

6 gll t21 t31 - 9 gll t21 t31 + 4 gll t31 + 2 gll t21 - 

2 2 3 3 3 

3 gll t21 - 8 gll t2l)/(t21 (t31 - t2l) ) 


s33 (t31 - g!3) (t31 - t21- gl3) (2 t31 - t21 - 2 gl3) 


2 2 
t31 (t31 - t21) 


sll (t21 + gll) (t31 + gll) (t31 + t21 + 2 gll) 

+ — — 

2 2 
t21 t31 
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s22 (t21 - gl2) (t31 + gll) (t31 - t21 + gl2 + gll) 


2 2 
t21 (t31 - t21) 

Do you want to continue displaying a ? 

6 

y; 


6 

2 2 

2 gll si (t21 + gll) (t31 + gll) (t31 + t21 t31 + gll t31 + t21 + gll t21) 


3 3 

t21 t31 

2 2 

+ 2 gll s2 (t21 + gll) (t31 + gll) (t31 - 2 t21 t31 + gll t31 + t21 

3 3 

- 2 gll t21)/(t21 (t31 - t21) ) - 2 gll s3 (t21 + gll) (t31 + gll) 


2 2 3 3 

(t31 - 2 t21 t31 - 2 gll t31 + t21 + gll t21)/(t31 (t31 - t21) ) 

2 2 2 2 
s33 (t31 - gl3) (t31 - t21 - gl3) s22 (t21 - gl2) (t31 + gll) 


2 2 2 2 
t31 (t31 - t2l) t21 (t31 - t21) 

2 2 

sll (t21 + gll) (t31 + gll) 


2 2 
t21 t31 
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Do you. want to display sll? 

y; 


y2 - 2 yl - 2 

sll = gll x2 (y2 - 1) y2 + gll xl (yl - 1) yl 

Do you want to display s22? 

y; 


y2 - 2 yl - 2 

s22 = gl2 x2 (y2 - 1) y2 + gl2 xl (yl - 1) yl 

Do you want to display s33? 

y; 


y2 

s33 = gl3 


2 yl 

x2 (y2 - 1) y2 + gl3 


2 

xl (yl - 1) yl 


(d6) done 
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APPENDIX F: Listing of Automatically Generated FORTRAN Code 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


This is the template subroutine to calculate tensor S and D. 
inputs are eigenvalues gll,gl2,gl3, and cmu(6) . emu is assumed to be 
engineering strain ( e) , e.g. the Cauchy-green deformation tensor 
cm(3,3) is related to emu (6) in the following fashion: 

cm(l , l)=cmu(l) , cm (2 , 2) =cmu(2) , cm(3,3)=cmu(3) , cmu(4)=2*cm(l ,2) , 
cm(5)=2*cm(2,3) ,cmu(6)=2*cm(l ,3) . 

The outputs are the second order tensor S(6) and forth order 
tensor D(6,6) are related in the following way: 

S=D*C 

S(1 , 1) - S(l) 

S(2,2) = S (2) 

S(3,3) = S (3) 

S(l,2) - S(4) 

S(2,3) = S (5) 

S (3 , 1) = S(6) 

C(1 , 1) = C(l) 

C(2,2) = C (2) 

C(3,3) = C(3) 

C(l,2) = C(4) 

C(2,3) = C(5) 

C(3 , 1) = C(6) 

subroutine compsd(gll ,gl2 ,gl3,cmu, s ,d) 
real*8 gll ,gl2,gl3,ts(3,3) ,td(3,3,3,3) 
real*8 delt(3,3) ,delt4(3,3,3,3) ,s(6) ,d(6,6) 
real*8 cmu(6) ,cm(3,3) 

converts cmu(6) to matrix cm(3,3) in a way that 
cm(l ,2)=cm(2 , l)=cmu(4) ,cm(2 ,3)=cm(3,2)=cmu(5) , 
cm(l,3)=cm(3,l)=cmu(6) 
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c 

do 5 i=i,3 
do 5 j=l,3 

if (i.eq.j) then 
iq=i 

cm(i, j)=cmu(iq) 
else if (i.ne.j) then 
if ((i+j) .eq.3) iq=4 
if ((i+j) . eq.4) iq=6 
if C (i+j) .eq.5) iq=5 
cm(i, j)=cmu(iq)/2 
end if 

5 continue 
c 

c Initiates the second identity tensor delt(3,3) which 

c is a 2X2 identity matrix 

c 

do 6 i=l,3 

delt (i ,i)=l .0 

6 continue 
c 

c Computes the forth order identity tensor delt4(3,3,3,3) 

c by definition, 

c 

do 7 i=l ,3 
do 7 j=l,3 

delt4(i , j , i , j ) =delt (i , i) *delt ( j , j ) +delt (i , j ) *delt ( j , i) 
delt4(i, j , j ,i)=delt4(i, j ,i,j) 

7 continue 
c 

c ****** *********************************************************** 
c For different eigenvalues gll t gl2,gl3 the computation is 

c different. 

c easel is gll#gl2#gl3 call subroutine comsdl. 

c case2 is gl3=gl2#gll or gll=gl3#gl2 or gll=gl2#gl3 then 

c call subroutine compsd2 . 

c case3 is gll=gl2=gl3 call subroutine ccmpsd3. 

c***************************************************************** 

c 
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if ((gll.ne.gl2) .and. (gl2.ne.gl3) .and. (gll.ne.gl3)) then 
call compsdMgll ,gl2,gl3,ae t,delt4,an,ts,td) 
else if (Cgl2.eq.gl3) .and. (gll .ne.gl3)) then 
call compsd2 (gll , gl2 , delt , delt4 , an, t s , td) 
else if ((gll.eq.gl2) .and. Cgl3.ne.gl2)) then 
gll=gl3 

call compsd2 (gll ,gl2, delt ,delt4,cm,ts,td) 
else if ((gll.eq.gl3) .and. (gl2.ne.gl3)) then 
gll=gl2 
gl2=gl3 

call coinpsd2 (gll ,gl2,delt ,delt4,cm,ts,td) 
else 

call compsd3 (gll , delt , del t4,ts,td) 
end if 
c 

c Rewrite the tensor ts(l,j) td(i,j,k,l) to S(i) and D(i,j) 

c respectively by using the symetric property, 

c converts ts(3,3) s(6) and td(3,3,3,3) to D(6,6) 

c 

do 8 i=l,3 
do 8 j=i,3 

if (i.eq.j) iq=i 
if (i.eq.l.and. j . eq.2) iq=4 
if (i. eq.2. and. j .eq.3) iq=5 
if (i . eq.l.and. j .eq.3) iq=6 
s(iq)=ts(i, j) 

8 continue 
c 

do 9 i=l,3 
do 9 j=i,3 

d(i,j)=td(i,i, j ,j) 

9 continue 
c 

do 10 i=l,3 

d(i ,4)=td(i ,i , 1 ,2)+td(i,i ,2 ,1) 
d(i,5)=td(i,i,2,3)+td(i,i,3,2) 
d(i,6)=td(i,i,3,l)+td(i,i,l,3) 

10 continue 
c 
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d(4,4)= (td(l,2,l,2)+td(i,2,2,l)+td(2,l,l,2)+td(2,l,2,l))/2 
d(4,5)= (td(l , 2,2,3) +td(l, 2,3,2) +td(2, 1,2,3) +td(2,l ,3,2))/2 
d(4,6)= (td(l ,2 ,l,3)+td(l ,2 ,3,l)+td(2,l ,l,3)+td(2 ,1 ,3, l))/2 
c 

d(5,5)=(td(2,3,2,3)+td(2,3,3,2)+td(3,2,2,3)+td(3,2,3,2))/2. 
d(5 ,6)=(td(2,3, 1 ,3)+td(2 ,3,3, l)+td(3,2,l ,3)+td(3,2,3,l) )/2 . 
d(6,6)=(td(3,l,l ,3)+td(3,l ,3,l)+td(l,3,l,3)+td(l ,3,3,l))/2. 
c 

do 11 i=l,6 
do 11 j=l ,6 

d(i , j) = d(j ,i) 

11 continue 

c 

c prints out the inputs gll ,cl2,gl3,cmu(6) 

c said outputs S and D 

c 

print*, * gll=’ , gll 

print*, * gl2= J , gl2 

print*, ’ glS^', gl3 

print*, * Input tensor 0(6):’ 

print*, (cmu(i), i=l,6) 

print*, "second order tensor S(6):" 

print*, (s(i), i=l,6) 

print*, "The forth order tensor D(6,6):" 
do 101 i=l ,6 
print*, (d(i,j) , j=l,6 
101 continue 

return 
end 


subroutine compsdl (gll , c312 ,gl3 , delt , delt4 , cm, ts , td) 
c 

real*8 gll ,gl2,gl3,ts(3,3) ,td(3,3,3,3) 
real*8 cm(3,3) ,delt(3,3) ,delt4(3,3,3,3) ,p(3,3,3,3) 
real*8 q(3,3,3,3) ,cnnn(3,3) , pi (3, 3, 3, 3) ,p21 (3,3 ,3 ,3) 
real*8 p3l(3,3,3,3) ,qll(3,3,3,3) ,ql2(3,3,3,3) ,p22(3,3,3,3) 
real*8 q2l(3,3,3,3) ,q22(3,3,3,3) ,a,b,c,al,a2,a3,a4,a5,a6 
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c 

c Obtains cmm(3,3)=cm(3,3)*cm(3,3) from subroutine product 
c 

call product ( cm ,anm) 
c 

c Uses the formula we derived in code to compute second order 
c tensor ts(3,3) . 
c 

do 25037 i=l ,3 
do 25038 j=l ,3 

ts(i , j)=a(gll ,gl2 ,gl3)*cmm(i, j)+b(gll,gl2 ,gl3)*cm(i , j) 
c(gll,gl2,gl3)*delt (i, j) 

25038 continue 
25037 continue 
c 

c call subroutine to compute the functions we defined 

c when we derived forth order tensor td, namely P(i,j,k,l) 

c and Q(l,j,k,l) which are the functions of cm(3,3) and 

c the matrix product 00111(3,3). 

c 

call pqcom ( cmm , cnm , pi , q) 
call pqcom(cmm,cm,p21 ,q) 
call pqcom(cm,cmm,p22,q) 
call pqcom(cm,cm,p31 ,q) 
call pqcom(cmm,delt ,p,qll) 
call pqcom(delt,cmm,p,ql2) 
call pqcom ( an, delt ,p,q21) 
call pqcom(delt,cm,P,q22) 
c 

c computes forth order tensor td(i,j,k,l) 

c 

do 25039 i-1,3 
do 25040 j-1,3 
do 25041 k=l ,3 
do 25042 1=1,3 

td(i,j , k , 1) =al (gll , gl2 ,gl3) *pl (i , j ,k,l)+ 
a2(gll ,gl2,gl3)*(p21(i, j ,k,l)+p22(i , j ,k,l))+ 
a4 (gll ,gl2 , gl3) *p31 (i , j ,k , 1) +a3(gll , gl2 , gl3) * 

(qll (i , j ,k,l)+ql2(i , j ,k,l))+a5(gll,gl2,gl3)* 

(q21 (i , j ,k,l)+q22(i , j ,k,l) )+a6(gll ,gl2 ,gl3) 
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25042 

25041 

25040 

25039 

c 

c 

c 

c 


c 


c 


*delt4(i, j ,k,l) 
continue 
continue 
continue 
continue 

return 

end 

a,b,c,al-a6 are the coefficients we derived in code. 

real*8 function a(gll ,gl2,gl3) 
real*8 gll ,gl2,gl3,sl ,s2 ,s3 

a=-sl (gll , gl2 ,gl3) / (-gll+gl2) / (-gll+gl3) +s2 (gll , gl2 ,gl3) / 
(-gll+gl2) / (-gl2+gl3) -s3 (gll , gl2 ,gl3) / (-gll+gl3) / (-gl2+gl3) 
return 
end 

real*8 function b(gll ,gl2,gl3) 
real*8 gll,gl2,gl3,sl,s2,s3 

b=s3 (gll ,gl2 , gl3) * (gll+gl2) / (-gll+gl3) / (-gl2+gl3) -s2 (gll , gl2 , 

gl3) * (gll+gl3) / (-gll+gl2) / ( - gl2+gl3) +sl (gll , gl2 , gl3) * (gl2+gl3) 

/(-gll+gl2)/(-gll+gl3) 

return 

end 

real*8 function c(gll,gl2,gl3) 
real*8 gll,gl2,gl3,sl J s2,s3 

c=-sl (gll , gl2 , gl3) *gl2*gl3/ (-gll+gl2) / (-gll+gl3) +gll* s2 (gll , 

gl2 ,gl3) *gl3/ (-gll+gl2) / (~gl2+gl3) -gll*s3(gll , gl2 ,gl3) *gl2/ 

(-gll+gl3)/(-gl2+gl3) 

return 

end 
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real*8 function ai(gll,gl2,gl3) 
real*8 gll ,gl2,gl3,sl ,s2 , s3, sll , s22 ,s33 

al=- sll(gll,gl2,gl3)/(-gll+gl2)**2/(-gll+gl3)**2+2.0* 
s2 (gll ,gl2 ,gl3) * (gll-2 . 0*gl2+gl3) / (-gll+gl2) **3/ 
(-gl2+gl3) **3-822 (gll , gl2 ,gl3) / (-gll+gl2) **2/ (-gl2+gl3) **2 
-s33(gll,gl2,gl3)/(-gll + g 13)**2/(-gl2+gl3)**2-2.0* 
si (gll , gl2 , gl3) * (-2 . 0*gll+gl2+gl3) / (-gll+gl2) **3/ 

( -gll+gl3) **3+2 . 0*s3 (gll , gl2 ,gl3) * (-gll-gl2+2 . 0*gl3) / 
(-gll+gl3) **3/ (-gl2+gl3) **3 
return 
end 

real*8 function a2(gll ,gl2 ,gl3) 
real*8 gll , gl2 , gl3 , si , s2 , s3 , sll , s22 , s33 

a2=s33 (gll , gl2 , gl3) * (gll+gl2) / (-gll+gl3) **2/ (-gl2+gl3) **2+ 

. s22 (gll ,gl2 ,gl3) * (gll+gl3) / (-gll+gl2) **2/(-gl2+gl3) **2+ 

. sll (gll ,gl2 ,gl3) * (gl2+gl3) / (-gll+gl2) **2/ (-gll+gl3) **2 
. -s3 (gll ,gl2 ,gl3) * (-2 . 0*gll **2-3 . 0*gll *gl2-2 . 0*gl2**2+ 

. 3 . 0*gll *gl3+3 . 0*gl2*gl3+gl3**2) / (-gll+gl3) **3/ (-gl2+gl3) **3 

. -s2 (gll ,gl2 ,gl3) * (2 . 0*gll**2-3 .0*gll*gl2-gl2**2+3 . 0*gll* 

. gl3-3 . 0*gl2*gl3+2 . 0*gl3**2) / (-gll+gl2) **3/ (-gl2+gl3) **3+ 

. sl(gll,gl2,gl3)*(-gll**2-3.0*gll*gl2+2.0*gl2**2-3.0*gll* 

. gl3+3 . 0*gl2*gl3+2 . 0*gl3**2) / (-gll+gl2) **3/ (-gll+gl3) **3 

return 
end 

real*8 function a3(gll ,gl2 ,gl3) 
real*8 gll , gl2 , gl3 ,81,82,83,311, s22 , s33 

a3=-sll (gll ,gl2 , gl3) *gl2*gl3/ (-gll+gl2)**2/ (-gll+gl3) **2-gll* 

. s22 (gll ,gl2 ,gl3) *gl3/ (-gll+gl2) **2/ (-gl2+gl3) **2-gll* 

. s33 (gll ,gl2 ,gl3) *gl2/ (-gll+gl3) **2/ (-gl2+gl3) **2+ 

s2 (gll ,gl2 ,gl3) * (gll**2*gl2-2 . 0*gll*gl2**2+gl2**3+gll**2 
*gl3-gll*gl2*gl3-2.0*gl2**2*gl3+gll*gl3**2+gl2*gl3**2)/ 

. (-gll+gl2) **3/ (- g 12+gl3) **3-sl (gll ,gl2 ,gl3) * (gll**3-2 . 0 
. *gll**2*gl2+gll*gl2**2-2 . 0*gll**2*gl3-gll*gl2*gl3+gl2**2*gl3 

. +gll*gl3**2+gl2*gl3**2) / (-gll+gl2) **3/ (-gll+gl3) **3-s3 (gll , 

. gl2,gl3)*(gll**2*gl2+gll*gl2**2+gll**2*gl3-gll*gl2*gl3+gl2**2* 

. gl3-2 . 0*gll*gl3**2-2 . 0*gl2*gl3**2+gl3**3) / (-gll+gl3) **3 

. /(-gl2+gl3)**3 
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real*8 function a4(gli ,gl2 ,gl3) 
real*8 gll ,gl2 ,gl3, si , s2 , s3 ,sll , s22 ,s33 

a4=-s33 (gll ,gl2 , gl3) * (gll+gl2) **2/(-gll+gl3) **2/ (-gl2+gl3) **2- 
s22 (gll ,gl2,gl3) * (gll+gl3) **2/ (-gll+gl2) **2/ (-gl2+gl3) **2 
-sll (gll ,gl2,gl3)* (gl2+gl3) **2/ 

(-gll+gl2) **2/ (-gll+gl3) **2+2 . 0*s2 (gll ,gl2 ,gl3) * (gll+gl3) * 

(gll**2-gll*gl2-gl2**2+gll*gl3-gl2*gl3+gl3**2) / (-gll+ 

gl2) **3/ (-gl2+gl3) **3-2 . 0*sl (gll »gl2,gl3) * (gl2+gl3) *(-gll**2- 

gll*gl2+gl2**2-gll*gl3+gl2+gl3+gl3**2)/(-gll+gl2)**3/(-gll+gl3 

) **3+2 . 0*s3 (gll , gl2 , gl3) * (gll+gl2) * (-gll**2-gll*gl2-gl2**2+gll* 

gl3+gl2*gl3+gl3**2) / (-gll+gl3) **3/ (-gl2+gl3) **3 

return 

end 

real*8 function a5(gll,gl2,gl3) 
real*8 gll , gl2 , gl3 , si , s2 , s3 , sll , s22 , s33 

a5=gll*s33 (gll , gl2 , gl3) *gl2* (gll+gl2) / (-gll+gl3) **2/ (-gl2+gl3) 

**2+gll*s22 (gll ,gl2 ,gl3) *gl3* (gll+gl3) / (-gll+gl2) **2/ (-gl2+gl3) 

**2+sll (gll , gl2 , gl3) *gl2*gl3* (gl2+gl3) /(-gll+gl2) **2/ (-gll+gl3) 

**2+s3(gll ,gl2,gl3)*(gll**3*gl2+gll**2*gl2**2+gll*gl2**3+gll** 

3*gl3+gll**2*gl2*gl3+gll*gl2**2*gl3+gl2**3*gl3-2.0*gll**2*gl3** 

2-5 . 0*gll*gl2*gl3**2-2 . 0*gl2**2*gl3**2+gll*gl3**3+gl2*gl3**3) / 

(-gll+gl3) **3/ (-gl2+gl3) **3-s2 (gll ,gl2 ,gl3) * (gll**3*gl2-2 . 0*gll 

**2*gl2**2+gll*gl2**3+gll**3*gl3+gll**2*gl2*gl3-5.0*gll*gl2**2* 

gl3+gl2**3*gl3+gll**2*gl3**2+gll*gl2*gl3**2-2.0*gl2**2*gl3**2+ 

gll*gl3**3+gl2*gl3**3)/(-gll+gl2)**3/(-gl2+gl3)**3+sl(gll,gl2, 

gl3) * (gll**3*gl2-2 . 0*gll**2*gl2**2+gll*gl2**3+gll**3*gl3-5 . 0* 

gll**2*gl2*gl3+gll*gl2**2*gl3+gl2**3*gl3-2.0*gll**2*gl3**2+gll 

+gl2*gl3**2+gl2**2*gl3**2+gll*gl3**3+gl2*gl3**3) / (-gll+gl2) **3/ 

(-gll+gl3)**3 

return 

end 
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real*8 function a6(gll,gl2,gl3) 
real*8 gll , gl2 , gX3 > si , s2 , s3 > sll » s22 , s33 

a6=-sll (gll ,gl2 , gl3) *gl2**2*gl3**2/ (-gll+gl2) **2/ (-gll+gl3) **2 
. -gll**2*s22 (gll ,gl2 ,gl3) *gl3**2/ (-gll+gl2) **2/ (-gl2+gl3) **2-gll 

. **2*s33 (gll ,gl2 ,gl3) *gl2**2/ (-gll+gl3) **2/ (-gl2+gl3) **2-2 . 0*gll 

. *s3 (gll , gl2 ,gl3) *gl2*gl3* (gll**2+gll*gl2+gl2**2-2 . 0*gll*gl3-2 . 0* 

. gl2*gl3+gl3**2) / (-gll+gl3) **3/ (-gl2+gl3) **3+2 . 0*gll *s2 (gll , gl2 , 

. gl3) *gl2*gl3* (gll**2-2 . 0*gll*gl2+gl2**2+gll*gl3-2 . 0*gl2*gl3+ 

. gl3**2) / (-gll+gl2) **3/ (-gl2+gl3) **3-2 . 0*gll*sl (gll ,gl2 ,gl3) *gl2* 

. gl3* (gll**2-2 . 0*gll *gl2+gl2**2-2 . 0*gll*gl3+gl2*gl3+gl3**2) / 

. (-gll+gl2)**3/(-gll+gl3)**3 

return 
end 

The sl,s2,s3,sll,s22,s33 are derivatives of W 

function si (gll ,gl2,gl3) 
common /param/ xl,x2,yl,y2 
real*8 xl,x2,yl,y2 
real *8 si ,gll ,gl2 ,gl3 

sl=2 . 0* (gll** (-1+yl) *xl*yl+gll** (-l+y2) *x2*y2) 

return 

end 

function s2(gll ,gl2,gl3) 
common /param/ xl , x2 , y 1 , y2 
real*8 xl,x2,yl,y2 
real*8 s2,gll,gl2,gl3 

s2=2 . 0* (gl2** (-1+yl) *xl*yl+gl2** (-l+y2)*x2*y2) 

return 

end 

function s3(gll ,gl2,gl3) 
common /param/ xl,x2,yl,y2 
real*8 xl,x2,yl,y2 
real*8 s3,gll ,gl2,gl3 

s3=2 . 0* (gl3** (-1+yl) *xl*yl+gl3** (-l+y2) *x2*y2) 

return 

end 
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function sll(gll ,gl2,gl3) 
common /param/ xl,x2,yl,y2 
real*8 xl,x2,yl,y2 
real*8 sil ,gll ,gl2 ,gl3 

sll=2 . 0* (gll** (-2+yl) *xl* (-1 . 0+yl) *yl+gll** (-2+y2) *x2* 

(-1.0+y2)*y2) 

return 

end 

function s22 (gll ,gl2 ,gl3) 
common /param/ xl , x2 , yl , y2 
real*8 xl,x2,yl,y2 
real*8 s22, gll,gl2,gl3 

s22=2 . 0* (gl2** (-2+yl) *xl* (-1 . 0+yl) *yl+gl2** (-2+y2) *x2* 

(-1 . 0+y2)*y2 

return 

end 

function s33(gll ,gl2,cl3) 
common /param/ xl,x2,yl,y2 
real*8 xl,x2,yl,y2 
real*8 s33,gll,gl2,gl3 

s33=2 . 0* (gl3**(-2+yl) *xl * (-1 . 0+yl) *yl+gl3** (-2+y2) *x2* 

(-1 .0+y2)*y2 

return 

end 

subroutine compsd2 (gll ,gl2 , delt , delt4 , cm, t s ,td) 
real*8 gll ,gl2,ts(3,3) ,td(3,3,3,3) 
real*8 cm(3,3) ,delt(3,3) ,delt4(3,3,3,3) ,pl(3,3,3,3) 
real*8 q2(3,3,3,3) ,ql(3,3,3,3) ,p(3,3,3,3) ,q(3,3,3,3) 
real*8 bl,b2,b3, abar.bbar 



c 


c Computes second order tensor ts(i,j) based on the formula 

c derived in code, 

c 

do 25043 i=l ,3 
do 25044 j-1,3 

ts(i,j)=abar(gll,gl2)*cm(i,j)+bbar(gll,gl2)*delt (i,j) 

25044 continue 
25043 continue 

c 

c Call subroutine to get P, Q which are defined in code, 

c 

call pqcom(cm,cm,pl ,q) 
call pqcom(cm,delt ,p,ql) 
call pqcom(delt ,cm,p,q2) 
c 

c Computes tensor td(i,j,k,l). 

c 

do 25045 1=1,3 
do 25046 j«l,3 
do 25047 k=l ,3 
do 25048 1=1,3 

td(i, j ,k,l)=bl(gll ,gl2)*pl(i, j ,k,l)+b2(gll ,gl2)* 

(ql (i , j ,k,l)+q2(i, j ,k,l))+b3(gll,gl2)*delt4(i, j ,k,l) 
25048 continue 

25047 continue 

25046 continue 

25045 continue 
return 
end 

c 

c abar ,bbar are bl, b2, b3 functions derived in code, 

c 

real*8 function abar (gll ,gl2) 
real8 gll,gl2,ssl,ss2 

abar=(-ssl(gll ,gl2)+ss2(gll,gl2))/ (-gll+gl2) 
return 
end 
c 
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real*8 function bbar(gll ,gl2) 
real*8 gll ,gl2,ssl ,ss2 

bbar= (-gll*ss2 (gll , gl2) +ssl (gll ,gl2) *gl2) / (-gll+gl2) 

return 

end 

real*8 function bl(gll,gl2) 
real*8 gll ,gl2, ssl ,ss2, sell ,ss22 

bl=2.0*ssl(gll,gl2)/(-gll+gl2)**3-2.0*ss2(gll,gl2)/(-gll+gl2) 
**3+ssll (gll , gl2) / (-gll+gl2) **2+ss22(gll ,gl2) / (-gll+gl2) **2 
return 
end 

real*8 function b2(gll,gl2) 
real*8 gll ,gl2,ssl , ss2 ,ssll ,ss22 

b2=-gll*ss22(gll,gl2)/(-gll+gl2)**2-ssll(gll,gl2)*gl2/ 

(-gll+gl2) **2-ssl (gll ,gl2) * (gll+gl2) / (-gll+gl2) **3+ss2 (gll , 

gl2) * (gll+gl2) / ( _ gll + gl2)**3 

return 

end 

real*8 function b3(gll,gl2) 
real*8 gll ,gl2,ssl ,ss2,ssll ,ss22 

b3=2.0*gll*ssl(gll > gl2)*gl2/(-gll+gl2)**3-2.0*gll*ss2(gll,gl2) 

*gl2/ (-gll+gl2) **3+gll**2*ss22 (gll ,gl2) /(-gll+gl2) **2+ssll (gll , 

gl2) *gl2**2/(-gll+gl2) **2 

return 

end 

ssl,ss2,ssll,ss22 are derivatives of W. 

function ssl(gll,gl2) 
common /param/ xl , x2 , y 1 , y2 
real*8 xl,x2,yl,y2 
real*8 ssl,gll,gl2 

ssl=2.0*(gll**(-l+yl)*xl*yl+gll**(-l+y2)*x2*y2) 

return 



function ss2(gll,gl2) 
common /param/ xl,x2,yl,y2 
real*8 xl,x2,yl,y2 
real ss2,gll,gl2 

ss2=2.0*(2 ,0*gl2**(-l+yl)*xl*yl+2.0*gl2**(-l+y2)*x2*y2) 
return 
end 
c 

function ssll(gll,gl2) 
common /param/ xl,x2,yl,y2 
real*8 xl,x2,yl,y2 
real*8 ssll,gll,gl2 

ssll=2 . 0* (gll** (-2+yi) *xl* (-1 . 0+yl) *yl+gll** (-2+y2) *x2* 
. (-1 .0+y2) 

return 
end 
c 

function ss22(gll,gl2) 
common /param/ xl , x2 , yl , y2 
real*8 xl,x2,yl,y2 
real*8 8s22,gll,gl2 

ss22=2 . 0* (2 . 0*gl2** (-2+yl) *xl* (-1 . 0+yl)*yl+2 . 0*gl2** 

. (-2+y2)*x2*(-1.0+y2)*y2) 

return 
end 
c 

subrout ine compsd3 (gl 1 , delt , delt 4 , t s , t d) 
real*8 gll,ts(3,3) ,td(3,3,3,3) ,delt(3,3) ,delt4(3,3,3,3) 
real*8 ccl,abbar 
c 

do 25049 i*l,3 
do 25050 j=l,3 

ts(i, j)=abbar (gll)*delt(i,j) 

25050 continue 
25049 continue 
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25054 

25053 

25052 

25051 

c 


c 


c 


c 


do 25051 i=l ,3 
do 25052 j-1,3 
do 25053 k=l ,3 
do 25054 1=1,3 

td(i , j ,k,l)=ccl(gll)*delt4(i, j ,k,l) 
continue 
continue 
continue 
continue 

return 

end 

real*8 function abbar(gll) 

real*8 gll 

abbar=sssl (gll) 

return 

end 

real*8 function ccl(gll) 

real*8 gll 

ccl=sssll(gll) 

return 

end 

function sssl(gll) 
common /param/ xl,x2,yl,y2 
real*8 xl,x2,yl,y2 
real*8 sssl.gll 

sssl=2 . 0* (3 . 0*gll** (-1+yl) *xl*yl+3 . 0*gll** (-l+y2) *x2*y2) 

return 

end 
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function sssll(gll) 
common /param/ xl,x2,yl,y2 
real*8 xl,x2,yl,y2 
real*8 sssll.gll 

sssl 1=2 . 0* (3 . 0*gll** (-2+yl) *xl* (-1 . 0+y 1) *yl+3 . 0*gll ** 

(-2+y2) *x2* (-1 . 0+y 2) *y2) 

return 

end 

This subroutine computes P and Q forth order tensors 
which we define in tensor D 

subrout ine pqcom ( cml , cm2 , p , q) 

real*8 cml (3,3) ,cm2(3,3) ,p(3,3,3,3) ,q(3,3,3,3) 
do 100 i=l ,3 
do 100 j=l ,3 
do 100 k=l ,3 
do 100 1=1,3 

p(i, j ,k,l)=cml(i,k)*an2(j ,l)+cml(i,l)*cm2(j ,k) 
q(i , j ,k,l)=p(i, j ,k,l)+cml(j ,l)*cm2(i,k)+cml (j ,k) 
*cm2(i,l) 

continue 

return 

end 

This subroutine computes matrix product cmXcm. 

subrout ine product (mat 1 , cmm) 
real*8 matl(3,3) ,cmm(3,3) ,sum 
do 25 i=l ,3 
do 25 j=l ,3 
sum=0 . 0 
do 26 k=l ,3 

sum= sum+mat 1 (i ,k) *matl (k, j) 
continue 
cmm(i, j)=sum 
continue 
return 
end 
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